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Abstract 

In this paper, we build on recent work using a mathematical programming ap¬ 
proach for incremental state update in analysis of non-linear mechanics models. 
In particular, we consider quasi-static analysis of continuum problems in the lin¬ 
earized kinematics regime, with non-linear material models described using con¬ 
vex energy functions. We find in this case that the classical displacement-based 
nested approach for incremental state update can be reformulated as solving a 
reduced dual optimization problem. This reformulation provides insights into 
the working of the algorithm, and eliminates the need for some heuristics. An 
important purpose of this paper is to further illustrate the unifying nature of 
the mathematical programming approach. We therefore present relationships 
with several of these types of algorithms recently presented in the literature for 
incremental state update. 

Keywords: mathematical programming; convex optimization; standard 
material; plasticity; Lagrange dual; state update 

1. Introduction 

Analysis of mechanics models with complex non-linear material behavior 
arises in different applications such as seismic response simulation of structures. 
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In a recent paper [1|, we discussed the manner in which when such material 


behavior is described using an energy approach [2l-l^ , state update of the model 
in each increment of a non-linear analysis can be cast into various mathemat¬ 
ical programs. Depending on the type of material model, the corresponding 
mathematical program could assume different forms such as convex optimiza¬ 
tion, complementarity problem etc. This mathematical programming approach 
represents an alternative to the classical displacement-based approach. In the 
latter, state update is carried out in a nested fashion with displacements com¬ 
puted at a global level, and stresses and other internal variables computed at 


a material-point level. It was also seen in [l| that models of non-linear be¬ 
havior that arose historically from disparate contexts, for example the Preisach 
model ^ and the Bouc-Wen model could be interpreted within the 

mathematical programming framework. 

In the present paper, we further explore the unifying nature of the mathe¬ 
matical programming approach. We consider quasi-static continuum problems 
in the linearized kinematics regime, with non-linear material behavior described 
by convex energy functions. In this case, we find that the classical displacement- 
based nested approach can be reformulated as solving a reduced dual optimiza¬ 
tion problem. Some differentiability results related to convex optimization prob¬ 
lems play an important role in developing this reformulation. This reformulation 
in turn provides insights into the working of the algorithm by means of some 
geometric constructions such as Figured and obviates some heuristics that are 
otherwise used. We envision that such insights will guide development of algo¬ 
rithms for more complex models, for example via successive convex programs 
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A particular case of non-linear material models with convex energy functions 
is rate-independent plasticity models. For this case, incremental state update 
at the material point level is known as return mapping, and has a well-known 
optimization format 15j. However, the optimization structure at the global 


level is not commonly recognized or utilized. The purpose of this paper is two¬ 
fold — (1) to reformulate the classical displacement-based nested approach for 
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incremental state update as an optimization problem (for quasi-static problems 
with convex energy functions), and (2) identify relationships with several other 
algorithms recently described in the literature. 

This paper is organized as follows. First in section [5J governing equations 
presented in [l| for single-degree-of-freedom systems with energy-based material 
models are extended to general spatially discretized systems. With some restric¬ 
tions, when discretized in time, these equations lead to an optimization problem 
in each time increments. In section |31 three forms of this optimization problem 
are discussed — a primal problem, a dual problem, and a reduced dual problem. 
The reduced dual problem is a reformulation of the classical displacement-based 
nested approach, and consists of solving an optimization sub-problem at the in¬ 
tegration point level. Solution of the reduced dual problem is the subject of 
section 21 Differentiability of the objective function is considered and deriva¬ 
tives are obtained, so that Newton’s method for unconstrained minimization 
can be applied as described in section 221 In section 21 some special classes of 
material models are discussed, where the integration point level optimization 
problem admits a simpler solution process. A numerical example is presented 
in section El to illustrate the working of the algorithm of section 14.21 Lastly in 
section [71 relationships between various solution strategies presented in recent 
years based on the mathematical programming approach are discussed. 


2. Governing equations and time discretization 


We begin by formulating the governing equations for a continuum model that 
has been spatially discretized, for example by the finite element element method. 
We employ a class of non-linear material models described completely by certain 


energy functions. As discussed in [l|, such a representation of material models is 


based on the generalized standard material framework [2,[l6|, and is also closely 


related to the hyperplasticity framework [3|, Il7| . In this energy framework, a 
material model is characterized by two convex (possibly nonsmooth) functions 
— a stored energy function ?/;(e, a), and a dissipation function (()(e, d), where e is 
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the strain, and a is a vector of kinematic internal state variables (such as equiv¬ 
alent plastic strain, damage etc.). We refer to the Fenchel-Legendre transforms 
of these functions as the complementary stored energy function '0'^(cr, C), and 
the complementary dissipation function where a is the stress and the 

( is the generalized stress conjugate to a. The governing equations discussed 
in [l| for simple single-degree-of-freedom dynamic models can be generalized to 
spatially discretized continuum models as 

momentum conservation : + V 


B^a 


= P 


element force 


inertia force 


deformation compatibility 


internal variable evolution 


diip^a, 0 


elastic deformation rate 


d2ij‘'{cr, C) 


reversible inter¬ 
nal variable rate 


damping force 

dlfp^'icr, C) 

^ ^ > 

plastic deforma¬ 
tion rate 

d2(p''{cr, C) 

V-,-^ 

irreversible 
internal variable 


Bv 9 0 


total deforma¬ 
tion rate 


9 0 


( 1 ) 

Here, v S K^dop jg vector of free velocity components, M, Cd € M^dopxaTdop 
are respectively positive semi-definite mass and damping matrices, and p € 
R'^DOP jg external load vector, a G is the stress vector comprised 

of Ncr components at each of the Nq material points. Similarly, ( G jg 

the vector of generalized stress internal variable. x R^c^g —R 

are the complementary stored energy and complementary dissipation functions. 
B G R^'^^gxVdop jg j;] 2 g Ijnearized deformation-displacement matrix. V denotes 
the gradient, and d denotes the subdifferential of a nonsmooth convex function 
[l^ . Vg and dq denote gradient and subdifferential with respect to argument 
q of the function. Equations ©2.3 are inclusions because '0'^ and (j)^ could 
be nonsmooth, and hence have set-valued derivatives. We make the following 
remarks on equation ©• 

1. As suggested by their format, equations o can be obtained as Euler- 
Lagrange equations of a generalized Hamilton’s principle 

2. The terminology we have used above in reference to spatial discretization 
follows the common setting, where displacements and velocities are asso- 
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dated with nodes, and stresses and other internal variables with element 
integration points. However, the formulation itself is not restricted to this 
type of spatial discretization. For example, in a stress-based finite element, 
cr would denote components of the stress approximation rather than the 
integration point values, or in an assumed-strain finite element, v would 
include components of the element-level assumed strain coefficients. 

3. For notational convenience, we take that equations (II|)2,3 have been multi¬ 
plied through by the appropriate integration weights and element Jacobian 
determinants. This explains our use of the term “deformation rate” in¬ 
stead of “strain rate” in equations 02.3, and the reason why the element 
nodal forces can be written simply as a in equation o 1 . Later in 
section ITTI the tangent stiffness matrix is of the form B^OB. 

4. In the above, we consider a restricted class of material models for which ip 

and (p are convex. One source of non-convexity is so-called non-associated 
flow in plasticity models. It can be shown that in some instances of 
non-associated flow models, the complementary dissipation function can 
be written as a) with explicit dependence on the kinematic inter¬ 


nal variables 


20(1 . Equations ([T]) then still apply, as discussed for simple 


single-degree-of-freedom systems in [l|. Explicit dependence on a then 
suggests solving the incremental state update as a succession of convex 
subproblems, computing a in an outer stage. 

5. When nonlinear kinematics is considered, the term B in equations © 1,2 is 
replaced by D£{u). £ is the nonlinear deformation-displacement map (i.e., 
the function that maps node displacements to element strains), u is the 
vector of free displacement components, and Dis the derivative operator. 
Thus B = D£’(0). In some instances, the kinetic energy in the first term 
of equation Oi may also be of the form M{u)v. This again suggests 

a successive convex programming approach, updating u in an outer stage. 

For the sake of concreteness, in what follows, we introduce the following addi¬ 
tional restrictions. 
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1. The complementary stored energy function ijj^ is smooth. Nonsmooth 
stored energy functions, such as those resulting from approximating tension- 
or compression-only behavior and some forms of damage, can however be 
treated in a similar manner to nonsmooth dissipation functions as de¬ 
scribed below [l|. 

2. The complementary dissipation function is of the form 

= Llc(o-,C) (2a) 

where C is a convex set given by 

C = {(ct,C)|</5(o-,C) < 0} (2b) 

Uc is its indicator function [l^ (Figure , and tp is a smooth vector-valued 
function, each component of which is convex. In plasticity models, (p is 
the yield function, and C is the elastic region. 

We note that this restricted framework is sufficient to describe nonlinear elasto- 
plasticity. We assume for simplicity of presentation that there is an identi¬ 
cal number Ny of yield functions at each material point. Thus (p : x 

]S,^(Ng With the restriction ([ 2 ]), the subdifferential of (/>'^ that ap¬ 

pears in equations © 2.3 can be written as (Figure fib)) 

C) = {A^V(^(cr, C)k(o-, C) < 0, A > 0, C) = 0} (3) 


For quasi-static analyses, the first two terms of equation o 1 are absent. 
Also using the fact that has been restricted to being smooth and equation 
©, equations o can be written in the alternate form 

B^a = p 
+ - Bv =0 

(4) 

^V2^^(a,C) + ATV2(^(a,C) =0 

C) > 0, A > 0, X^picr, 0 = 0 


6 




CC 



(a) Indicator function of a convex subset 
C of 



(b) Subdifferential of the indicator function of C; The subdifferential is shown 
at three points — xi in the interior if C, X 2 where the boundary is smooth, 
and X 3 where the boundary is not smooth 


Figure 1: Indicator function and its subdifferential for a convex set C = {x G 

R^\(p^ix)<0,i = l,2} 


Equations (U) may be formally discretized in time. Using Backward Euler 
with a time increment At, we get 






,C+i)-ViV’^(a",D 

At 

At 


+ (A"+i) 
+ (V+i) 


> 0, A”+^>0, 


7 / 1 _ q ,. 


= 0 
= 0 
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where u S K^dop jg vector of free displacements. Multiplying equations 
©2,3 by At, rearranging, and introducing the symbols /r ;= — it" and 

A := At A, equations (O become 

ViV'"(a"+\ C”+^) + A'^Vi(^(a"+\ C”+^) - {B^i + C")) = 0 

^ 

V2i/>"(a”+\C+^) + A^V2(^(CT”+\C”+^)- V2V'"(a",0 =0 

'-^ 

-■be. 

= p”+i 

yj(a"+\C”+^) > 0, A>0, A^yj(a"+\ C"+^) = 0 

Equations ([HI) become the starting point for optimization problems of different 
forms as discussed in the next section. We note that prescribed displacements 
can be accounted for by simply replacing all occurrences of the term (Bfi + ba-) 
above by {Bji + ^prsc^prsc _|_ where is the increment of prescribed 
displacements, and B^^^^ is the deformation-displacement displacement matrix 
associated with DOF with prescribed displacement. 


3. Optimization problems 

3.1. Primal optimization problem 


We recognize immediately that equations 
optimality (KKT) conditions 


are the Karush-Kuhn-Tucker 


2l| for the following minimization problem. 


) = argmin ^ C) - 
subject to B^a = 

< 0 

where ba and b^ are as defined in equation ([6]). The incremental displacements 
are the Lagrange multipliers corresponding to the equilibrium constraints (mo¬ 
tivating the notation p). The Lagrange multipliers A corresponding to the yield 
constraints can be interpreted as incremental equivalent plastic strains. 
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o is a convex minimization problem, since the objective function and in¬ 
equality constraints are convex, and the equality constraint is linear [^|. We 
refer to this as the primal optimization problem. We note that corresponds 


to the mathematical program referred to as “approach 1” in 
to equation (46) of reference 


Q. 


It is identical 


23| (where the notation k is used for C), and sim¬ 


ilar to equation (5) in 2^. It is also the same as equation (17) in [ 23 , where 
dynamics is considered as well (and the notation {F,() is used for (cr, C))- We 
assume in the following that 

1. © has a solution. 

2 . tjj^ is strongly convex, so that the solution of © is unique, and the min 
in equation (HU is meaningful. 

3. © satisfies Slater’s constraint qualification condition [22|, i.e. that it has 
a strictly feasible point. 


We write © in a more convenient form by noting that the objective func¬ 
tion and the constraints are separable over integrations points. The objective 
function can be written as 

Ng 

n (^J C) ^ ^ V’ (f^miCm) ~ ^am'^rn ~ (8) 

m—1 

and the optimization problem as 


(^„+l^^„+l) ^ ^(( 7,0 

Ng 

subject to ^ B^Um = P (PRIMAL) 

m—1 

‘fi(crTn,Cm) <0, TO = 1, . . . , A^G 

where subscript to denotes components of a vector or rows of a matrix corre¬ 
sponding to integration point to. We have used and ip in the above equa¬ 
tions instead of and pm, to minimize notational clutter. Models where the 
material is not homogenous, so that these functions are different at different 
integration points presents no additional difficulty. We have also dropped the 
superscript n -|- 1 on the load vector p for brevity. 
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The Lagrangian associated with (IPRIMAL|) is 

/ Ng \ Ng 

-^((c) ^ (a)) = I X! - p| + ^m<P(o-m,Cm) 

\ m=l j m—1 

Ng 

= // p ^ ^ [^0 {o'rm Cm) H“ '^mV^(^m 5 Cm) H“ ^(rm) ^^mCm] 

m=l 

A^g 

= + X! ( Cm ) ’ ( aD 

m=l 

(9) 

where we have dehned 

(( Cm ) ’ ( Am )) ■” V” ('^mi Cm) T Cm) {Bmf^ + b^rm) O'm ^CmCm 

(10) 

The Karush-Kuhn-Tucker (KKT) optimality conditions for (IPRIMALD follow 
from the Lagrangian ([9]) as 

^ = 0 : Di '(/'‘"(o'm, Cm) + Dl ip{am, Cm) “ + bamV =0, m = 1, . . . , TVg 

^ = 0 : D 2 Ip''{am, Cm) + D 2 ip{am, Cm)-bJ^ = 0 , m = 1 ,..., Nq 

Ng 

BmUm = P 

m—1 

Cm) ^ 0; ^m ^ 0; '^mV^(*^m; Cm) — 0; ^ — 1; ■ • ■ ; -^G 

(KKT) 

These are merely equations ([ 6 ]) separated over the integration points. Since 
the gradient is the transpose of the derivative, while equations ([ 6 ]) 1^2 are in 
terms of column vectors, equations (IKKTD i 9 are in terms of row vectors. We 
next construct the dual of problem (IPRIMALI) , which in section |4] leads to a 
reformulation of the classical displacement-based approach as an optimization 
problem. 

3.2. Dual optimization problem 

Duality is familiar in mechanics from the classical principles of total potential 
energy and total complementary potential energy, and from the upper- and 
lower-bound theorems of limit analysis. Associated with every optimization 
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problem is a dual problem, which often has a useful interpretation. For convex 
problems such as (IPRIMAL|) that satisfy constraint qualification conditions as 
assumed above, strong duality holds. This means that a dual optimal solution 
is equal in objective function value to a primal optimal solution, and that the 
primal-dual solution pair satisfies the KKT conditions. The KKT conditions 
are also sufficient conditions for a solution of the primal problem 221. 


To construct the dual problem for (IPRIMALI) . we start from the Lagrange 


dual function 


22|, defined as the minimum of the Lagrangian of 11'^ over (cr,(). 


:=mm£((^) ,(^)) ( 11 ) 

(<^,0 

We are able to write min here instead of inf, because of the strong convexity 
assumption on Let = 1,...,Nq be the minimizers in (ED. 

Functions X) and (^(fj.,X} are then implicitly defined by the following 

optimality conditions. 

A), C(m, A)) + A^ ^(a;;(/i. A), C(/i, A)) - = 0 

Lb A), C(p, A)) + AT D, A), C(m, A)) - =0 

( 12 ) 

i.e., given (/r, A), the functions cr^(/J-, A) and A) can be evaluated by solving 
equation ED- Clearly, this is equivalent to unconstrained minimization of the 
Lagrangian function of equation ED- 

The Lagrange dual function 11 of equation (ITT]) may be written explicitly as 

Ng 

n(/i,A) = -/rTp- ^ [V’"«(^,A),C(M,A))+AT^(a;;(A.,A),C(/^,A)) 

m—1 

+ berm) A) — b^enCmif^^ A) 

(13) 

The dual optimization problem corresponding to dPRIMAL]) is then 


minn(/i. A) 
subject to A > 0 


(DUAL) 


As discussed above, solving this problem is completely equivalent to solving 
(|PRIMALI) . (IDUALI) is a convex optimization problem with simple bound con- 
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straints. It is closely related to the mathematical program referred to as ap¬ 
proach 2 in [l|. In section next, we reduce this to an unconstrained op¬ 
timization problem in ^ alone. This forms the basis for reformulation of the 
classical displacement-based approach in section H) 


3.3. Reduced dual problem 

Due to the separability of the objective function 11, minimization over A in 
(IDUALI) mav be carried into the summation in m- Thus defining the function 


Ng 


n(/i) := -pJp- ^ 


mm 

A,n>0 L 


A), C(/i, A)) + A), C(/i, A)) 

~{Bml^ 3-berm) ^) ~ ^(mCmil^T ^) 

(|DUAL|) can be reduced to minimization of /r alone as 


(14) 


minn(^) (REDUCED) 

This reduced problem leads to the reformulation of classical displacement-based 
approach as an optimization problem. We note that in this way, the solution 
of (IDTTAIi is carried out in a nested fashion — the minimization over Am > 0 
within the summation is at the integration point level, while the minimiza¬ 
tion over n is at the global level. We reiterate that the integration point-level 
minimization format is well-known in the context of so-called return-mapping 
schemes, discussed further in section [S] However, the minimization structure at 
the global level (IREDUCEDI) is not commonly recognized. We take on the solu¬ 
tion of (IREDUCEDI) in section |4l In the remainder of this section, we develop 
an explicit expression for H. 

Consider the minimization within the summation in equation (1141) . 

min A), A)) -b A^v?(cr;;(Ai, X),Cmil^, A)) 

- (Bmfi + b.m)^a*e^it,, A) - bJ^Cmil^, A) (IPDUAL) 

subject to Am > 0 
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where the prehx IP is used in the equation tag since the minimization is at 
the integration point level. Let Am(/r) be the minimizer. We recognize that 
(|IPDUALI) is the Lagrange dual of 


min C) - (-Bm/i + barnVa - 


(IPPRIMAL) 


subject to (p{a, C) < 0 


Let (dm(/r), Cm(M)) be the minimizer, given implicitly by the optimality condi¬ 
tions 


Di Ip {cFyyi (/i), (^)) “t“ \rn (A^) IA (/^) ; (m)) ^(rm) — 0 

lA (o'm(M), ^ , CmifJ')) “ = 0 


= 0 


(15) 


By strong duality, the optimum values of (IIPDUALI) and (IIPPRIMAL|) are 
identical. Therefore, LI can be written as 



m—1 


(16) 


We note as an aside that solve the KKT conditions of 

(|IPPRIMAL|) . and that the identities XmifJ-)) = ^mip) and Cm(M: XmifJ-)) = 
hold. Equation (fTBl) provides an explicit expression for LI. In the next 
section, we present an algorithm for its minimization, i.e. to solve (IREDUCEDI) . 

4. Optimization reformulation of the classical displacement-based nested 
approach 

In this section, we consider numerical solution of the reduced dual problem 
(|REDUCEDI) . We seek to use Newton’s method for unconstrained optimization. 

Eor this, we require the first and second derivatives of the objective function (ITBl) . 

We discuss differentiability of this function, and computation of these derivatives 
in section Q] below. We then present Newton’s method in section Solving 
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the reduced dual problem in this manner may be interpreted as a reformulation 
of the classical displacement-based nested approach as an optimization problem. 
The reason for this interpretation is also discussed in section [4.II 

4.I. Differentiability and derivatives 

We seek to obtain the first and second derivatives of If needed for mini¬ 
mization using Newton’s method. In computing these derivatives, is helpful to 
think of (lIPPRIMALI) as a parametric convex optimization problem with pa¬ 
rameter p. The first derivative of n is well-defined, and can be obtained by 
differentiating m as 

A^g 

= Y, - P (17) 

m—1 

The optimality condition for (IREDUCEDI) . Vli(/i) = 0, is thus 

A^g 

YB^am{p)=P (18) 

m—1 

which is the equilibrium equation. 

Equation (fTSl) can also be obtained directly from (IS|), since m 1 ^ 2,4 are the 
KKT conditions of (lIPPRIMALI) . When derived this way however, the mini¬ 
mization structure (IREDUCEDI) at the global level is not obvious. Recognizing 
this structure enables direct application of Newton’s method for unconstrained 
minimization in section EM without need for heuristics, often used for exam¬ 
ple in step-length determination. The classical displacement-based nested ap¬ 
proach used (1181) as the starting point, rather than the minimization problem 
(IREDUCEDI) . It is in this sense that we consider (IREDUCEDI) as a reformula¬ 
tion of the classical approach. 

To compute the second derivative of LI, we need the derivative of the func¬ 
tions dm(,p)- However, in trying to differentiate dTSl) to obtain this, we realize 
that drn{p) and Cmip) are not differentiable when they are degenerate solutions 
of (lIPPRIMALI) . i.e., when ip{am{p)Xm{p)) = 0 and \m{p) = cQ. In fact at 

^Indeed these degenerate points feature in the context of the first derivative (|17ll as well, 


14 

















such points, these functions are only directionally differentiable Q . When the 
minimizer of (IREDUCED|) is such a degenerate point, the rate of convergence 
of Newton’s method is only super-linear rather than quadratic 271. 


Here, we proceed formally ignoring these degenerate points, the set of which 
has measure zero in . We first introduce the index set 


T — 

-^m. — 


{ke{i,...,Ny}\iXM)k>o} (19) 

The optimality conditions (IlSp can then be rewritten as 

Dl 5 Cm (m)) + m (/r))fc Di (fikid- m ifk) 5 Cm (/i)) “h b(jrn^ — 0 


k£Xrr 


1^2'4^ (5'm(/^)i Cm(/^)) H" E(^ m Cm(/^)) 

k^Xm 

= 0, for k Glm 

Differentiating this, we get 


= 0 


K.. 



KJ, 

KCC 




0 




1 V 


D(Am(/r)x„:) 



Bm 

= 

0 


. 0 , 


(20a) 


where 


K(T(T — ^ ( (Atti (/i))A; E^ (dm (/r), Cm (/r)) 

feGim 

) Cm (m)) + E 0 m {n))k E>2 Di (pk{d- m ifk) ) Cm (m)) 
kGXm 

l^CC = 1^(dm (m), Cm (A^)) + E Om(Ai))fcE^‘A’fc(dm(Ai),Cm(M)) 

(DLV5(dm(Ai),Cm(M))), 

(^i:^(A’(dm(Ai),Cm(Ai))), 


feGX„ 




In the last two equations above, we have used MATLAB-style indexing 


28|. 


The subscript Im ■ extracts the submatrix consisting of rows Im and all the 


where an argument has to be made that Am(/r)^ D(^ = 0 at such points. 
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columns. 

The derivative can be obtained by solving (I20al) . It is of the form 

KmBm, where is referred to as the consistent tangent stiffness matrix 291. 


When the submatrices in (I20al) have a simple form, it may also be less expensive 
to compute the consistent tangent directly as 


= K- 


K,, + K„^K7}kJ. - $ 


where 


XCC = Krr - KJ^RzIK, 


^T^-l 




KZ 


‘■CC 


crC^aa^<y< 


5 = 


(21a) 


(21b) 


Now, by differentiating equation m, the second derivative of the objective 
function 11, the tangent stiffness matrix, may be written as 


Ng 


V2n(/r) = ^ 


( 22 ) 


m—1 


4-2. Optimization algorithm 

With the first derivative Vn(^), the vector of unbalanced forces, and second 
derivative V^n(/r) the tangent stiffness matrix, at hand, the energy n(/i) can be 
minimized using the standard Newton’s method for unconstrained minimization 
[^,2, outlined in Procedure [TJ A Newton direction Sp is computed by solving 


the linear system 


= -Vn(/r“"^) 


(23) 


To assure convergence for large increments, a line search is typically needed. 
A step length s is computed by backtracking to satisfy the following sufficient 
decrease condition. 


(24) 


n(/r + sSp) < n(^) + /3svn(/i)^(5^ 

where /3 is a sufficient decrease parameter to be in the interval (0,0.5) 
in the example in sectionlHlas 10“'^. Solving the problem as one of optimization 


2l|, taken 
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removes the need to use any other heuristics to ensure convergence. The working 
of the backtracking procedure can be visualized as shown in Figured in section 

E 


Procedure 1 Newton’s method for minn(/r) 

1: Given force p, prescribed displacement /rP''®'’, {ba-m, b^rn) for each integration 

point m, and starting point 
2: for iter ^ 0, MAXITER do 

3: Compute and using Procedure^ 

4: if ||Vn(/j,’’-^®’’’)|| < tol then 

5: break 

6: end if 

7: Compute search direction 5^ by solving (|23)) 

8: s = 1 

9: for nback MAXBACKTRACK do \> backtracking line search 

10: Compute + s5y) using Procedure [2] 

11 : if sufficient decrease (l24ll is obtained then 

12: break 

13: end if 

14: s i— "fs [> 7 is a backtracking parameter 

15: end for 

16: ^ + s6n 

17: end for 


5. Solution of (jIPPRIMALIl for specific constitutive models 

In the framework described here, a constitutive model is specified by a com¬ 
plementary stored energy function '0®, and a yield function tp or complemen¬ 
tary dissipation function 0®. In general, (lIPPRIMALI) may be solved (step H] 
of Procedure ED using any standard algorithm for inequality constrained con¬ 
vex optimization such as an interior-point method or a sequential quadratic 
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Procedure 2 Evaluate n(/i), Vn(^) and V^n(/i) 
1: Given /i, p, 

2: n ^—pjp ; vn ^—p ; v^n ^ o 

3: for m ^ 1, Nq do 

4: Compute {am, Cm) by solving (lIPPRIMALIl 


5: 

6 : 

7: 

8 : 


n- 


{am. Cm) {^m^ “f ^crm) am ^QmC'i 


Assemble B{^am into VII 


Computer Km using either of equations (I20al) or (I21al) 


Assemble BmKmBm into V^II 


9: end for 


programming method 221. Alternately, it may also be advantageous on oc¬ 
casion to solve (IIPDUALII . which has simple bound constraints. Such a general 
approach may often be necessary, particularly for more complex models such 
as multi-surface plasticity models. However in some cases, specific forms of the 
complementary stored energy function and the yield function ip may facil¬ 
itate simpler strategies to solve (lIPPRIMALIl . These strategies are known as 
return-mapping algorithms. In the following, we discuss one such case, namely 
an elastoplastic model characterized by isotropic linear elasticity, linear kine¬ 
matic hardening, nonlinear isotropic hardening, and von Mises yield condition. 
This model is described by 


C) = V -f ^C^hH ^Ckh + V'ih(Cih) 

C) = Y/('^“Ckh)^P(u-Ckh) - {ay + Cih) 


(25) 


Here, Ckh and Cih are components of the internal variable C corresponding to 
kinematic hardening and isotropic hardening respectively. Ckh is often referred 
to as the back stress, and is the same type of object as a. Ckh is a scalar. C and H 
are matrices of elastic and kinematic hardening moduli. The quadratic terms in 
'0'^ correspond to linear elasticity and linear kinematic hardening. The function 
is not necessarily quadratic, and represents nonlinear isotropic hardening. 
The relationship between this form of 0'^ and the hardening model discussed in 
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Q is discussed in Q . von Mises yield function in both 3D and 2D plane-stress 
cases can be expressed in the form tp above with appropriate choice of the matrix 
P (Tabled]) [l5|. ay is the uniaxial yield stress. 

The optimality conditions (ITSl) corresponding to the model (ESI) are 



(26) 


V’icr, C) < 0, A > 0, X(p{a, 0 = 0 


where the prime on denotes the first derivative. For brevity, we have written 
a for dm{p) etc. and dropped the subscript m. Since there is a single yield 
condition, there are only two possibilities. 

1. (/3((T, 0 ^ 0. In this case, the material point is elastic in the increment, 
and A = 0 from equation ([26ll /i. From equations (IMl) i 9 y. an elastic trial 
state is computed. 



(27) 


If < 0 then this is the solution of (lIPPRIMALIl . 

2. If (^(cr*’''®'\ ^triai^ ^ concluded the material point plastifies. So 

A > 0, and (Eni)i,2.3, 4 constitute four equations in the four unknowns cr, Ckh, 
Cih and A. To solve these, a further property of the matrices involved may 
be invoked. In both 3D and 2D plane-stress cases, the matrices C, H and 
P share the same eigenvectors. Therefore, there is an orthogonal matrix 
Q (Tabled])) consisting of these eigenvectors as columns, which simulta¬ 
neously diagonalizes all three matrices, i.e. PQ is diagonal matrix Ap 
etc. This has been pointed out for the 2D case in 0 . The 3D case is 
generally presented differently, by decomposing the stress into volumetric 
and deviatoric parts. In both cases, however, Q represents transformation 
of stress components to volumentric-deviatoric coordinates. Equations 
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Table 1: von Mises yield condition matrix P and diagonalizing matrix Q 


1 , 2 , 3,4 can be rewritten as 
A 


Ckh — 


3 i^y + Cih) 
A 


-CP{a - Ckh) - = 0 


3 

C f 


■ Cih) 


i?P(fT-Ckh)-crh'" =0 


(28) 


V'fh(Cih)-A^--5c,,=0 

\/(ci - Ckh)^^’(CT - Ckh) = \/f ('^y + Cih) 

where is the component corresponding to Cih- Introducing the 
coordinate transformation a = Qua- and Ckh = Qy subtracting (l28l) 9 
from (I28l) i . and multiplying through from the left by , we have 

\ 

A 


/ 


2/o’ VChh ~ 


V 


3 (^y P Cih) 


^ic+H)P Q ' - CS"0 (29) 


where the diagonal matrix A(^c+h)p = (C + H)PQ. Solving equation 
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(1^ ^ for A and substituting, we get 


^ (Cih) ^Cih 

2 CTy + (^ih 


-1 


y<T - 2/Ckh =1^+2 


^{C+H)P 


iT/^trial /-trial\ 

-Ckh ) 


Lastly, substituting this in (|281 ) a and squaring, we get the following scalar 
equation for Cih. 



(30) 


where r goes from 1 to 6 in 3D and from 1 to 3 in 2D. Equation (15(1)) can 
also be written equivalently as 



(31) 


where Oih is the kinematic internal state for isotropic hardening, and ijAh 
is the Legendre transform of 'i/'fh- o^ih is the equivalent plastic strain [l|. 
Equation m is a scalar non-linear equation in Oih, which can be solved by 
Newton’s method. The other states can then be calculated as summarized 
in Procedure m 

The consistent tangent can be calculated using m- The above process of 
solving (lIPPRIMALI) for the constitutive model (1251) is summarized in Procedure 


El 


We recount that a simpler strategy to solve (lIPPRIMALI) for the constitutive 
model (|25l) is enabled by the following features. 

1. The model has only one yield condition 

2. a and C are uncoupled in the complementary stored energy function '0'” 

3. The part of C that appears in a non-quadratic manner in namely 
the isotropic hardening state Cm, is a scalar. This helps in reducing the 
problem to a scalar equation as in dSU. 
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4. C, H and P share the same eigenvectors, so that there is a diagonalizing 
coordinate transformation Q. 


In addition, in the 3D case, (Ap)i = 0 since the von Mises yield constraint 
applies to the deviatoric part of the stress, and {Ai^c+h)p) 2:6 = 2(G + Gkh) 
where G and Gkh are the elastic and kinematic hardening shear moduli. The 
solution of (lIPPRIMALI) can also be interpreted as a radial projection of the 
trial state on the elastic region, and is referred to as radial return map ping 
[l5 . In the 2D plane stress case, squaring the yield condition as done in ^ 
further simplifies the expression for the consistent tangent. However, here we 
keep the yield function as shown in (1251) 9. so that it remains a convex function 
and so that the 3D and plane stress cases can be treated in a uniform manner. 
We also note that if isotropic hardening were also linear, then the objective 
function in (lIPPRIMALI) can be written as 


1 


(a-cr*™yG-i(a-. 


'*) + ^(C - 


Thus the solution can be interpreted as the closest point projection in the 
norm of the trial state on the elastic region. In the absence of such 
simplifying features, (jlPPRIMALI) must be solved by a general constrained con¬ 
vex optimization strategy. 


6. Numerical example 


In this section, we present a numerical example to illustrate the application of 
ProcedureHJ We use a finite element model of a plate with a circular hole subject 
to uniform extension at the end, considered by Simo and Taylor 30|. Due to 


symmetry, only a quarter of the plate needs to be represented as shown in Figure 
[21 In order to compare the numerical solution and convergence characteristics, 
we use the same displacement-based finite element and discretization as used in 
this reference. We note however, that the formulation presented in this paper 
can be used with other finite element types. A constitutive model characterized 
by linear elasticity, linear isotropic hardening, and a single von Mises’ yield 
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Procedure 3 Solution of (IIPPRIMALI) for constitutive model ((^ 
(return map) 

2 ; Compute = cr" + CBn, and = Cm 

3; if (/?(cr‘™',C*™0 < 0 then 
4: cr = cr*"*^', C = A = 0 

5: K = C 

6 : else 

7 : Compute QT (^trial _ ^trial) 

8 : Solve (IHTl) for Oih 

9- Cih — l/^ih (rTih) 

10 : Calculate A using (l28l) -i 

11 : Compute cr and Ckh from equation ( 12 ^ 1 -? 

12 : Calculate K using equation (EU 

13: end if 


function, is used to represent material behavior. This is a special case of the 
material model (|25]) with V'fh(Cih) = 5^~^Cih- Thus the Procedure [3] can be 
used to solve (IIPPRIMALI) . Material nrooerties are chosen followinff H as 
E = 70, v = 0.2, cTy = 0.243 and h = 2.24, where E is Young’s modulus and 
V is Possion’s ratio. Extension is applied to the plate by imposing a uniform 
displacement at the top edg e of the plate. Of the three sequences of increments 


considered in reference 


30(1 . we only present results for the most challenging 


case, where two displacements increments of size 0.5 are applied in succession. 
In Procedure[2l the backtracking parameters are taken as /3 = 10“^ and 7 = 0.5. 


Table [5] shows norms of the residuals that arise in applying the optimization 
algorithm of Procedure [T] for each increment. The approximate doubling of 
significant digits in the residuals in later iterations suggests quadratic rate of 
convergence. The table also shows the step length used in each iteration, from 
which we see that at most two backtracking steps used and that the full Newton 
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Figure 2: Finite element model of a plate with circular hole 

step is always used in the last several iterations. The performance is comparable 
to that reported in [^, with fewer backtracking steps used attributable to use 
of the sufficient decrease condition. The von Mises stress computed and the 
spread of plasticity are shown in Figure |31 

A typical backtracking process used to obtain the step length is depicted in 
FigurelH Such plots represent cross-sections along the Newton search path, and 
show the following as functions of the step length s. 

1. Objective function, II(/x -|- s6fi) 

2. First order approximation, n(/i) -I- sVn(/r)^(5/r 

3. Second order approximation, li(/.t) -t- 6 fj, + 

4. Sufficient descent criterion, n(^) -|- /3sVn(^)^(5^ 

These plots illustrate the working of the backtracking search, as well as serve to 
check derivative computations during implementation. In Figure |4al the iterate 
is far away from the solution. Thus the objective function deviates considerably 
from its second order approximation, and a step length of 0.5 is required to 
satisfy the sufficient decrease condition (l24l) . On the other hand in Figure 
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Increment (Step Length) 


Iteration 

I 


2 


I 

6.66e+00 


5.93e+00 


2 

2.I8e+00 


1.44e+00 

(0.5) 

3 

1.38e+00 

(0.25) 

6.80e-01 

(0.5) 

4 

1.57e+00 

(0.5) 

7.00e-01 

(0.5) 

5 

1.33e+00 

(0.25) 

5.33e-01 

(0.25) 

6 

l.I2e+00 

(0.25) 

3.43e-01 

(0.5) 

7 

8.59e-0I 


1.34e-01 


8 

3.96e-0I 

(0.5) 

3.I2e-02 


9 

I.75e-0I 


9.63e-03 


10 

2.06e-02 


8.87e-04 


II 

I.53e-03 


2.45e-06 


12 

6.15e-06 


1.89e-ll 


13 

I.42e-I0 





Table 2: Norm of of unbalanced force, ||Vn(^)|| 





(a) 


(b) 


Figure 3: Numerical example (a) von Mises stress, V Pa, after increment 
2 (b) Spread of plasticity (green — yielded in increment 1; red — yielded in 
increment 2; orange — some integration points yielded in increment 2) 
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[ibl the iterate is close the solution, the second order approximation closely 
follows the function, and the full Newton step is taken. We note that because 
of the optimization reformulation, we are able to select the step length based 
on sufficient decrease of the objective function, rather than on heuristics such 
as number of iterations. In this example, despite the rather large displacement 
increments, at most two backtracking steps are needed in any iteration as seen 
from Tabled] 


7. Relationships with other approaches 


In section[3l we derived three equivalent optimization problems for incremen¬ 
tal state update — a primal problem dPRIMAL]) . the dual problem (IDUALI) . 
and a reduced dual problem (IREDUCEDI) . In the present section, we briefly 
review relationships these problems share with other approaches recently pre¬ 
sented in the literature. 

There has been recent work on solving the primal problem in the context of 
quasi-static elastoplastic models. Krabbenhoft et al. 23| arrive at the primal 


problem through a different route, and solve it directly by an interior point 
method. Bilotta et al. 2^ also solve the primal problem directly, but by a 

solves a KKT 


sequential quadratic programming (SQP) approach. Wieners H 
system similar to ([S)) by SQP. When is quadratic (corresponding to linear 
elastic and linear hardening models), and when the yield function is of the form 
of a 2-norm (as is the case with the von Mises function in section Ej), then 0 
can be recast as a second order cone program (SOCP) This can in turn 


be reformulated as a semi-definite program (SDP) [33|, which can be solved by 
efficient interior point methods that have been recently developed. Krabbenhoft 
et al. Q,Q also use this approach to approximate a non-associated flow model 


by a convex one. Sivaselvan et al. 


251 solve a primal-type problem that arises 


dyni 

i3. 


in dynamics problems. Further relationships to dynamic problems are discussed 
In all of these primal approaches, a concept similar to the Sherman- 


Morrison-Woodbury formula is used to reduce the linear system to be solved 
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(b) Increment 2, iteration 10 


Figure 4: Visualization of backtracking for step length determination 


25l| . Bilotta et al. 2^ discuss the relationship between the primal optimization 
problem and the displacement-based nested approach through the dual problem, 
albeit lightly differently than done in this paper. In the KKT conditions 


are solved directly as a mixed complementarity problem (MCP). 

As noted in section ITTl equation (fT^ can also be obtained directly from 
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Using 

projection 

[31] 



Nested approach 
by locally computing 
plastic strains [39] 


Figure 5: Interrelation of various approaches to elastoplastic problems using 
optimization 


since (l 6 |)i, 2,4 are the KKT conditions of (lIPPRIMALI) . However, when viewed 
in this manner, the minimization structure at the global level is not apparent. 
In the context of plasticity models, , Cmi^J-)) is referred to as a return 

map. As discussed at the end of section [SJ for certain classes of problems, the 
return map defined by (lIPPRIMALI) can be viewed as a projection [l^, l^ |. 
Other approaches are also possible, for example local problem is considered as 
a linear complementarity problem (LCP) in 


^ Following the seminal paper 
Newton’s method has been successfully applied to solve (TT^ in many cases. 
However, since the return map is not differentiable, convergence properties of 
Newton’s method cannot be rigorously established. Consequently, semismooth 
methods have been presented together with convergence analyses j^ . 

Optimization problems for quasi-static elastoplasticity can also be derived 
starting from variational inequality formulations presented for example by Han 
and Reddy 3^ . An optimization problem developed in this manner in terms of 


displacements p and kinematic internal states or plastic strains a is considered 
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by Alberty et al. 


. They use the optimization structure to develop error 


401 de¬ 


measures and an adaptive mesh-refining strategy. Gruber and Valdman 
rive a minimization problem of the form (IREDUCED|) for a class of elastoplastic 
models again starting from a variational inequality. They prove the differentia¬ 
bility of a function similar to If using the Moreau-Yosida theorem. Recognizing 
that the second derivative is not continuous, as discussed in section im they use 
a slant Newton method. Some of the interrelations described above are shown 
schematically in Eigure[SJ 


8. Concluding remarks 

We take a mathematical programming approach to the incremental state up¬ 
date of non-linear mechanics models with material behavior described by stored 
energy and dissipation function. In the case where these energy functions are 
convex and of a specific form (equation ([2|)), the classical displacement-based 
nested approach can be reformulated as a reduced dual optimization problem 
(|REDUCEDI) . We have presented an optimization reformulation of the classical 
displacement-based nested approach to incremental state update. This reformu¬ 
lation allows visualizing the working of the algorithm by means of geometrical 
constructions such as FigurelH and further illustrates the unifying nature of the 
mathematical programming approach to state update. Connections with related 
algorithms recently presented in the literature have also been discussed. We en¬ 
vision that such understanding gained will aid in developing algorithms for more 
complex non-linear material models, such as those involving non-associated flow, 
softening etc., for example by means of successive optimization strategies. 


NOTATION 

B Linearized strain-displacement matrix (g 

Bm Rows of the matrix B corresponding to integration point m 
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Deformation-displacement matrix associated with DOF with 
prescribed displacement 

Cd Damping matrix (s R^dopxWdop) 

C Matrix of elastic moduli in the constitutive model equation 

(|^ 


C Convex elastic region 

D, Dq Derivative of a differentiable function; subscript q denotes with 

respect to argument q 

E Young’s modulus (used in the numerical examples) 

£ Nonlinear deformation-displacement map used only in the re¬ 

marks following equation 0 

G, Gkh Elastic and kinematic hardening shear moduli used only at the 

end of section [5] 

H Matrix of kinematic hardening moduli in the constitutive model 

equation (USD 

I Identity matrix of appropriate size 

2m Index set of activated yield conditions defined in equation (I19|) 

IFcrcr, ^CC Derivative matrices dehned in equation (I20b|) 

Km Consistent tangent stiffness matrix for integration point m 

K(^i^ Derivative matrix defined in equation (I21bl) 

C Lagrangian of the primal optimization problem (equation @) 

Cm Contribution to primal Lagrangian of integration point m (equa¬ 

tion (fTUl) ! 

M Mass matrix (g R-^DOPxAfDOP) 

MAXBACKTRACK Maximum number of backtracking steps in Procedure [1] 
MAXITER Maximum number of Newton iterations in Procedure [T] 

-/Vdof Number of free degrees of freedom 

Nq Number of integration points 
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Number of stress components per material point (for example, 
3 in 2D problems and 6 in 3D problems) 

Number of internal variables per material point (for example, 4 
for 2D elastoplasticity with one yield condition and combined 
isotropic and kinematic hardening) 

Matrix in von Mises yield condition (051) 9 
Derivative matrices defined in equation (I21b|) 

Diagonalizing coordiirate transformation matrix (Table [T]) 
Time increment 

Diagonal matrix containing the eigenvalues of the symmetric 
matrix □, used in equation (OH)) 

Derivative matrices defined in equation (I20bl) 

Objective function of dual problem, defined in equation (ED, 
explicit formula in equation (IT3l) 

Objective function of the reduced dual problem defined in 
equation (ED 

Objective function of the primal optimization problem (equa¬ 
tion dH)) 

Derivative matrix defined in equation (I21bl) 

Kinematic internal variable at a material point 
Isotropic hardeing internal state conjugate to (jkh, used in equa¬ 
tion (IHT|) 

Sufficient decrease parameter used in choice of step length 
(equation (ED) 

Newton search direction (equation (1231) ') 

Strain at a material point 

Backtracking parameter used in Procedure [1] 

Vector of equivalent plastic strain increments for the entire 
model (s or for a single material point (£ R^y) 








A, A”+i 


Vector in equations m and ® (e corresponds to equiv¬ 

alent plastic strain rate; superscript denotes time increment 
index 

Am Minimizer of (lIPDUALj) 

/i Incremental displacement, — it" 

^ in iteration number iter of Newton’s method (Procedure [T]) 
Increment of prescribed displacement 
v Poisson’s ratio (used in the numerical examples) 

(j) Dissipation function for a material point 

(jf Complementary dissipation function for a material point, or 

for entire model, depending on context 
ip Yield function 

tp Stored energy function for a material point 

Complementary stored energy function for a material point, or 
for entire model, depending on context 

tjj^^ Part of complementary stored energy representing kinematic 

hardening in the constitutive model equation (I25I) 

V'ih Legendre transform of used in equation (I3I|) 

cr, cr", CTm Stress at a material point (g or collection over all ma¬ 
terial points (g depending on context; superscripts 

denote time increment index, and subscripts material point in¬ 
dex 

dm Minimizer of (lIPPRIMALI) 

Elastic trial stress (equation (071) 1 

tTy Uniaxial yield stress in von Mises yield condition (051) 9 

am Minimizer of equation m defined implicity in equations m 

C) Cm Generalized stress internal variable at a material point (g R^*^ ) 

or collection over all material points (g depending on 
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Ckh; Cih 

context; superscripts denote time increment index, and sub¬ 
scripts material point index 

Components of internal variable C corresponding to kinematic 

and isotropic hardening, used in constitutive model equation 

(|^ 

c 

Sm 

Minimizer of equation dill) defined implicity in equations dl2|) 


(m Minimizer of (lIPPRIMALI) 

^triai^ ^tnai^ ^triai giastic trial State (equation ([?7|1 

ba-m Defined in equation ([5]); components of ba- corresponding to 



integration point m 

Defined in equation ([S]); components of bi^ corresponding to 

integration point m 

^Ckh; ^Cih 

Components of 5,^ corresponding to Ckh and Cih, used following 

equation (05)) in section [5] 

h 

Isotropic hardening modulus used in the numerical example 

iter 

Iteration count in Newton’s method in Procedure [1] 

k 

Yield function component index within an integration point 

m 

Integration point index G {I,..., Nq} 

n 

Time increment index 

nback 

Number of backtracking steps in Procedure [1] 


External load vector, and its values at time n -I- 1 (G 

r 

Index used in equation (1501) 

s 

Length of Newton step in[I] 

t 

Time 


tol Convergence tolerances for Newton’s method in Procedure [T] 

u, u", Displacement at free DOF, and its values at times n and n + 1 

(e 
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Velocities at free DOF, and its values at times n and n + 1 

(e 

X 

Generic variable used in Figure [1] 

Vfyi VCkh 

Gomponents of stress and kinematic hardening internal state 

in volumetric-deviatoric coordinates 

<1 

<1 

Gradient of a differentiable function; subscript q denotes with 

respect to argument q 

d, dg 

Subdifferential of a nonsmooth convex function; subscript q 

denotes with respect to argument q 

Uc 

Indicator function of the convex elastic region C 
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